The Cosmology of f(R) Gravity in the Metric Variational Approach 
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We consider the cosmologies that arise in a subclass of f(R) gravity with f(R) = R+fi 2n+2 / (-R) n 
and n £ (—1,0) in the metric (as opposed to the Palatini) variational approach to deriving the 
gravitational field equations. The calculations of the isotropic and homogeneous cosmological models 
are undertaken in the Jordan frame and at both the background and the perturbation levels. For the 
former, we also discuss the connection to the Einstein frame in which the extra degree of freedom 
in the theory is associated with a scalar field sharing some of the properties of a 'chameleon' field. 
For the latter, we derive the cosmological perturbation equations in general theories of f(R) gravity 
in covariant form and implement them numerically to calculate the cosmic-microwave-background 
temperature and matter-power spectra of the cosmological model. The CMB power is shown to 
reduce at low Vs, and the matter power spectrum is almost scale-independent at small scales, thus 
having a similar shape to that in standard general relativity. These are in stark contrast with what 
was found in the Palatini f(R) gravity, where the CMB power is largely amplified at low Vs and the 
matter spectrum is strongly scale-dependent at small scales. These features make the present model 
more adaptable than that arising from the Palatini f(R) field equations, and none of the data on 
background evolution, CMB power spectrum, or matter power spectrum currently rule it out. 
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I. INTRODUCTION 

Theories of f(R) gravity as an explanation of the dark 
energy problem have attracted much recent attention. 
These theories are of particular interest since f(R) mod- 
ifications to general relativity (GR) appear naturally in 
the low-energy effective actions of the quantum grav- 
ity and the quantization of fields in curved spacetimc. 
These theories are also conformally related to GR with 
a self-interacting scalar field |l[ and both the early time 
inflation and the late-time acceleration of the universe 
could be resulted by a single mechanism in such theo- 
ries. In refs. 0, Q > the specific model in which the cor- 
rection is a polynomial of R 2 , R ab R a b and R abcd R a bcd is 
considered and the analysis there shows that late-time 
accelerating attractor solutions exist. Meanwhile, mod- 
els with R, R ab R a i, corrections are discussed within the 
Palatini approach 0, H, H, 0], in which the field equa- 
tions are second order, and similar acceleration solutions 
are found. The Palatini-/(i?) theory of gravitation was 
then tested using various cosmological data such as Su- 
pernovae (SN), Cosmic Microwave Background (CMB) 
shift parameter, baryon oscillation and Big Bang Nucle- 
osynthesis (BBN), in d S, EE E3, E3- Recently, con- 
straints from CMB and matter power spectra have also 
been given in [l3|, EH EH • These constraints (especially 
that from the matter power spectrum) successfully ex- 
clude most of the parameter space, making the model 
indistinguishable from the standard ACDM in practice. 
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Turning to the metric f(R) gravity theories, the field 
equations become fourth order and more difficult to han- 
dle [l6j |. Until now, much attention has been focused on 
the solar system tests of the theory and the existing re- 
sults appear to exclude a viable f(R) cosmological model 
that leads to the current cosmic accelerating expansion 
(see, for example i _El El, El, US HH HI and references 
therein; also see (23| and references therein for arguments 
against these). Regardless of these local considerations, 
however, we can look at the cosmological behavior of the 
theory in order to arrive at an independent test. The 
cosmological constraints on f(R) gravity is relevant also 
because one can simply imagine that the baryons do not 
see the modifications to GR (a possibility suggested in 
[HI) evades the solar system tests in a somewhat artifi- 
cial way), or transform to the Einstein frame and con- 
sider a scalar field model mathematically equivalent to 
the Jordan frame f(R) gravity [26j], but with the scalar 
field coupling only to the dark matter. For further dis- 
cussion of the cosmology of the f(R) dark energy models 
see, for example, [13, HH, H^] and references therein. 

In a recent study, [24J, the authors find that the class 
of Friedmann cosmological models arising in a theory 
with f(R) = R + ^ n+1 ^/{-R) n and n > (note the 
difference from their original paper because of our sign 
convention, R < 0), which was thought to lead to the 
late-time acceleration because the correction term to GR 
becomes significant at late times, cannot reproduce the 
matter-dominated era of conventional cosmology (yield- 
ing a oc t 1 / 2 rather than a oc t 2 ^ 3 where a is the scale 
factor and t the cosmic time) and so will be ruled out by 
measurements of the redshift distance (see also [IBJ). In 
this work they do not consider the case of — 1 < n < 
which is able to give a standard matter-dominated era, as 
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we show below. In a later work, [30|, these authors con- 
sider the cosmological viabilities of general f(R) models 
and have included this possibility. 

Since there are few known models with simple forms 
of f{R) having the capability to consistently describing 
the whole evolutionary history of the universe (see how- 
ever [3l| for a discussion on this), we will also consider 
the perturbation evolutions of the model. To this end 
we derive the covariant and gauge invariant perturba- 
tion equations and put them into the public CAMB code 
[32| in order to obtain the CMB and matter power spec- 
tra. Note that the dynamics of perturbations for f(R) 
gravity have also been considered in [H, HI] . In [33| the 
background expansion is fixed to match ACDM while in 
[3^ ] two specified models are considered which do not 
give standard matter-dominated eras and the calculation 
is done in the Einstein frame; thus both analyses differ 
from the work reported here. 

This paper is organised as the following: in Section HT1 
we briefly introduce the cosmological model and list the 
background and perturbation evolution equations that 
are needed in the numerical calculation. In Section IIIII 
we solve for the background evolution of the model nu- 
merically, and incorporate the perturbations equations 
into CAMB in order to calculate the CMB and matter 
power spectra. The discussion and conclusions are then 
presented in Section HVl Throughout this work we will 
set c = h = 1 and only consider the case of a spatially flat 
universe filled with photons, baryons, cold dark matter 
(CDM) and three species of effectively massless neutri- 



II. FIELD EQUATIONS IN F(R) GRAVITY 

In this section we very briefly summarise the main in- 
gredients of f(R) gravity in metric approach, and list the 
general perturbation equations that govern the dynamics 
of small inhomogcncitics in this theory. 



A. The Generalised Einstein Equations 

The starting point for the metric-/(i?) gravity is the 
Einstein-Hubert action, 



S = 



in 



(1) 



in which k — 8ttG with G being the gravitational con- 
stant and R = R(g a b) is the Ricci scalar. Varying the 
action with respect to the metric g ab gives the modified 
Einstein equations 

FR ab - \gabf{R) + (Sa&V c V c - V Q V h )F - nT ab , (2) 

where F = F(R) = df(R)/dR and T ab is the energy- 
momentum tensor. The trace of Eq. ([2]) reads 



with T° = T = p — 3p. Let us call this the structural 
equation, which relates R (or F) directly to the energy 
components in the universe. However, the term V c Vc-F 
enters here which makes the equation a dynamical rather 
than algebraic one. Recall, in contrast, that in the Pala- 
tini f(R) gravity the metric g ab and connection T^ c are 
treated as independent variables with R = R(g ab ,T^ c ); 
the variation of the action is taken with respect to both 
of these two variables with the resulting field equations 
being second order and the structural equation simply 
an algebraic relation. Now we can also make a conformal 
transformation, 



gab 



gab, 



(4) 



from the original Jordan frame to the Einstein frame (the 
tilded quantities are the Einstein frame ones from here 
on), in which we obtain the following action for the metric 
f(R) gravity: 



S = 



2k 



- V(tp) + 



(5) 



In Eq.02 we have defined a new scalar field ip = V6uj/k = 
v /3/2KlnFand its potential V = (f-FR)/2nF 2 . Quan- 
tities in the Jordan and Einstein frames are related to 



each other by [24 1 
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p = pe ,p 



pe 



-4u> 



, dt = e u dt, a — e"a. 



(G) 



We see that in the Einstein frame the scalar field couples 
minimally to gravity but couples conformally to matter, 
and that the generality in the functional choice of f{R) 
manifests itself in the generality of the potential V (</?). 

Since we implement our measurements in the Jordan 
frame, we shall treat the Jordan frame metric g ab as the 
physical one. Hence the difference between /(i?) gravity 
and GR (f(R) — R) can be understood as a change in the 
way that the spacetime curvature, and thus the physical 
Ricci tensor R a b, responds to the distribution of matter. 



B. The Perturbation Equations 

The perturbation equations in general theories of f(R) 
gravity derived in this section which uses the method of 
3+1 decomposition [33, Hfl- For more details, we refer 
the reader to these references and only briefly outline its 
main ingredients before listing out results. 

The main idea of 3 + 1 decomposition is to make space- 
time splits of physical quantities with respect to the 4- 
velocity u a of an observer. The projection tensor h ab is 
defined as h ab = g ab — u a ii b which can be used to obtain 
covariant tensors perpendicular to u. For example, the 
covariant spatial derivative V of a tensor field Tj.'.'.'" is 
defined as 



FR - 2/ + 3V C V CJ F = kT, 



(3) 



v«rj;;. 



hfK 



hZS7*Tt 



(7) 
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The energy-momentum tensor and covariant derivative 
of 4- velocity are decomposed respectively as 

Tab = TTafc + 2<7( a Ufc) + pU a U b - ph ab , (8) 

V a u 6 = <r ab + m ab + ^9h ab + u a A b . (9) 

In the above, n ab is the projected symmetric trace-free 
(PSTF) anisotropic stress, q a the vector heat flux vec- 
tor, p the isotropic pressure, a ab the PSTF shear tensor, 
^ab = VfaUw, the vorticity, 9 — V c w c = 3a/ a (a is the 
mean scale factor) the expansion scalar, and A b = u b 
the acceleration; the overdot denotes time derivative ex- 
pressed as cj) = u a V a <t>, brackets mean antisymmetrisa- 
tion, and parentheses symmetrization. The normaliza- 
tion is chosen as « a u„ = 1. 

Decomposing the Riemann tensor and making use of 
the modified Einstein equations, Eq. @ , we obtain, after 
linearization, five constraint equations 



3F - _ 



3F 2 







tab + ^ 


2F J 






B ab + | 







4F 2 

Sab - V C S d ( a £ fe)ec d 



<^ab 



0;(20) 



C _ d e 
'd(a s b)ec u 



-^^ C ^d( a e b)t 



0,(21) 



where the angle bracket means taking the trace-free part 
of a quantity. 

Besides the above equations, it is useful to express the 
projected Ricci scalar R into the hypersurfaces orthogo- 
nal to u a as 



R 



k(p + 3p) - f 





ng a 
F 

B ab 



V b £ afc 



V»B ab 



V c {e a \ d u d w ab ); 
9V a F V a F 2V a 9 
3F F 3 
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1 • •• V 2 F 

-(F^ + 3F) + - r (22) 



(10) 

V b a a6 + V b ti7 ah (ll) 



The spatial derivative of the projected Ricci scalar, 7] a 
iaV a i?, is then given as 



t b)ec u > 



(12) 



1 










V b TT ab + 


(I- 




2F K 







2 - 



a ~ a 



1 



VaF 



2a 



3F 
2F 



a ■ a - - n 

-eV a F - -V a V 2 F, 



V Q 
(23) 



1 

2F 2 



K (p + p)^aF FS7 a F + FV a F (13) anc j j tg propagation equation by 



1 

— > 

2F 



V c g d + (p + p)vj cd 



cd ni b 



(14) 



Here, e a bcd is the covariant permutation tensor, £ ab and 
S a f, are respectively the electric and magnetic parts of 
the Weyl tensor W abc di given respectively through £ ab — 



29 a - - 9 



F 
F 



aV Q V • A 



u c u d W ac bd and 6 ah = 



?it it e. 



2"* "* °ac ' *efbd- 

In addition, we obtain seven propagation equations: 



p+ (p + p)0 + V% 

So + g^«a + (P + P)^a - V a p + V b 7T af) 



0;(15) 
0;(16) 



-!^V a V-g-|;V a (V 2 F)-. (24) 

As we are considering a spatially flat universe, the spa- 
tial curvature must vanish on large scales which means 
that R = 0. Thus, from Eq. ([22]). we obtain 
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3_F 







V 2 F np 
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2 

a afc + - 




3F 
" 4F 


^ab ~ 




2F K7rab " 





9 - V a A a 



J_ 

2F 



V< Q A b) 



«(P + 3p) 



3 



10- 



/ ^ 
— + — f 
F F 



-F 



0. (25) 



0;(17) 



This is just one of the modified Friedmann equations of 
the metric f(R) gravity, and the other modified back- 
ground equations (the other Friedmann equation and the 
energy-conservation equation) can be obtained by taking 
the zero-order parts of Eqs. (TT71 fT5|) . as 



3F 



«P ,J_ 
F 2F 



m+-9w-S7 [a A b] = 0;(19) 



T^ab 



0; (26) 
P+(p+p)9 = 0. (27) 



1 








3F\ 


K 

2F 


K a b + 


(I- 


2F j 


1 

~2F 


K 




h P)<Jab 


+ V 



It is easy to check that when f(R) = R, we have F = 1 
and these equations reduce to those of GR. 

To conclude this section, we want to point out another 
way to derive the above perturbation equations. It fol- 
lows by treating the new contributions in the Einstein 
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equation from the f(R) gravity modification as an ad- 
ditional effective energy momentum tensor [3||. In this 
way one can write 

Gab = K T a b 

= ^T a b + jV a V b F~^(f + K T-V 2 F)g ab , 

where an overbar represents the total effective quantity, 
total effective energy momentum tensor in this case. Now 
using the relations 

p = T ab u a u\ 

V = -\h ah T ab , 

q a = h d a u c T cd , 
ir ab = h c a h%T cd + ph ab 

we can identify the components of the total effective en- 
ergy momentum tensor as 



The subsequent results follow exactly those well known in 
GR, just with the components of the energy-momentum 
tensor being replaced by the effective ones. We have 
checked this approach leads to the same set of perturba- 
tion equations as we gave above. 



III. NUMERICAL RESULTS 

This section is devoted to the numerical calculations we 
have performed. First, we obtain the background evolu- 
tion of the present model, which we discuss both in the 
Einstein frame and in the Jordan frame, and compare 
the two. Once the background evolution at some time is 
known, we can implement our set of perturbation equa- 
tions on the CAMB code, with the background quantities 
at arbitrary times (as required by the code) obtained by 
interpolations. Some of the details are presented below. 



A. The Chameleon Mechanism in f(R) Model 

The starting point for obtaining the background evo- 
lution is the Eqs. (|25l - [27)) . However, before looking at 
this calculation, we discuss the evolution of background 
quantities in the Einstein frame, which provide us with 
a clearer physical picture. 



From the variation of the action, Eq. ([5]), we can obtain 
the following field equations: 



(28) 



<p + 3H<p + V v = \hiP-Zp), 



p + 3H(p + p) 



-{p-3p)tp, (29) 



where H = da/adt, quantities p and p include contribu- 
tions from both matter and radiation components, and 
V^p = dV/dip. For the model under consideration, f(R) 
is given by 



f(R) = R 



/'■' 



2(n+l) 



(-R)n 



so we have 



2 \ "+ 1 



F(R) 



Kp — 


2F ( ^ 


f 3F 1 
^-{f + 2F + 2F^ F ~ 


\- F6), 






up = 


K 

2F {p - 


^ + JF-^ 2F + F9) - 


F 

2F 1 


and 




nq a = 


1 


-lv a F-±6V a F, 






V{<p) 


KTT ab = 


1 

-pK7T a b 


+ j V (o V 6>J F + -a ab . 










1 + n 



exp 



(30) 



7, ex P ~^V6ntp 

2k V 3 



exp 



(4V 



1 



(31) 



Here, our definitions ensure that [exp(v / 6K<p/3) — l]/n 
is always positive. At early times, when p?/(—R) <C 1, 
we have F — > 1, <p — > 0; thus it follows immediately that 
when — f < n < we have, in this limit, F — > f ~, tp — > 
0~, [exp(V6K<p/3) — l]/n — > + and so V — > oo. On the 
other hand, if tp — > — oo, we also have V — > oo. Also, the 
potential has a global minimum, 



V min (ip) = 



JL In 2 "+ 2 
2k n+2 ' 



/' 



8/4(71 + 1) 



(n + 2) 



i+2 



at <pa - 

Since there is a coupling to the matter, the evolution 
of tp is not determined solely by its potential, but by an 
effective potential Ves(tp) given by (from Eq. (|28|)) 



, . ... , 1 / 2V6^ x 
V e ${tp) = V(p) + -p m exp I —p 



in which p m is the energy density of the non-relativistic 
matter components, where we have used the Jordan- 
frame matter density, which is independent of p, and 
the fact p m — 0,p ra d = p r ad/3- The new minimum of the 
effective potential, tpo, satisfies p>Q < po < 0. From ip 
rightwards, the effective potential is dominated by V{tp), 
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which is essentially a runaway one in this region; from 
ipo leftwards, it is dominated by the matter coupling, 
which increases exponentially as ip becomes more nega- 
tive. Thus, this situation has the advantageous proper- 
ties of a Chameleon field [13, l38l. |39l. lib! ] , the cosmology 
of which has been studied in |4ll. 142]. However, there 
is one important difference between our case and theirs. 
In order to see this, we show that the effective mass- 
squared of small oscillations around the potential mini- 
mum, m 2 = d 2 V eS {if )/dip 2 , is given by 



np n 



( Kp n 



n+1 



v 3n(n + l) V M 2 

when v^M is small enough (note that n < 0). Mean- 
while, with v^l^l <C 1 the quantities in the Jordan and 
Einstein frames are essentially equal and we have 



""ip 
H 2 



( Kp n 



n(n + 1) V M 2 



n+l 



(32) 



in which Q m = np m /3H 2 . During the matter-dominated 
era f2 m ~ 1 and in the radiation dominated era O m cx a 
while (Kp nl /fi 2 ) n+1 cx a~ 3 (™ +1 ). Thus, for n > -2/3, the 
ratio m^/H 2 always increases with increasing redshift, 
and, because {np m / p 2 ) n+1 ~ 1 today, it is much larger 
than 1 at early times, just as in the standard Chameleon 
scenario. The differences lies in the fact that, at late 
times, in the present model m^/H 2 ~ 1 for |n| not too 
close to 0. 

We can now depict the evolution history of the field 
as follows: at the early times the effective mass of <p is 
much heavier than H so that the field is attracted to- 
wards its effective potential minimum, oscillates around 
it, but with the amplitude the oscillations being gradu- 
ally damped so that it eventually tracks the minimum. 
A similar analysis to the one in ref. [4l| can be made to 
show that the field rolls slowly along this attractor. As 
the universe evolves, m v / H gets smaller so that eventu- 
ally the field begins to lag behind its effective potential 
minimum and behaves like a normal quintessence field, 
which contributes a dynamical dark energy. In the far 
future, when p m — ► 0, the dynamics of the field is de- 
termined by its potential only, so that it would evolve 
towards the potential minimum, ipo, and stay there, af- 
ter which the universe commences de Sitter expansion in 
both the Einstein and Jordan frames. 

At this stage, we can also look at what happens if 
n = 0~. Obviously, from Eq. fl5D]) , F — 1~ and thus 
the field ip is confined to 0~ for all of the cosmic history. 
From Eq. (j3"2"|) . m 2 /i? 2 — * oo so that any deviation from 
- decays immediately. The potential V(<p) is, however, 
nonzero according to Eq. (|31|) ; in fact, in this limit we 
have 



lim V 



2k' 



So, this represents a non-dynamical field with constant 
potential, which is essentially a cosmological constant, 



and in this case we recover the correct ACDM limit, as 
expected. 

Finally a comment on the solar system constraints on 
the f(R) models. As discussed by the authors of Ref. [49| . 
when the baryons are allowed to see the f(R) modifica- 
tion to GR, the parameter region where one has a dynam- 
ical dark energy does not overlap with that in which the 
chameleon effect could successfully suppress the scalar- 
tensor type deviation from GR in the solar system [Zj|. 
Here, as mentioned in Sec. HI however, we are merely con- 
cerned with employing a "parametrized post-Friedman" 
description [33j of the cosmological effects of the f(R) 
model and assume that these deviations from the stan- 
dard gravitational model do not affect the baryons. Fur- 
thermore, as noticed in [l9j . the exclusion of f(R) mod- 
els using solar system constraints relies on the assump- 
tion that the transition from GR-dominatcd dynamics to 
scalar-tensor dynamics on these scales occurs adiabati- 
cally, which has not been investigated in detail. The non- 
uniqueness of the static spherically symmetric solutions 
to higher-order gravity theories is also a complicating fac- 
tor when determining the solar system bounds because 
of the absence of a Birkhoff theorem; see Refs. [3 EH 
for more details. 



B. Background Evolution 

Although we rely on the transformation to Einstein 
frame to obtain a physical picture for the present model, 
our numerical calculations for the background evolution 
will be implemented in the Jordan frame to be consistent 
with later perturbation calculations. We describe these 
calculations in this section. 

Following [3(|, we make the following definitions 



x 3 



F 
~~FH' 
f 

6FH 2 ' 
R H 



X4 — f^rad — 
RFr 



6H 2 H 2 



+ 2, 



3FH 2 ' 



F ' 
-RF _ X3 

/ x 2 ' 



(33) 
(34) 

(35) 
(36) 
(37) 
(38) 



where F R = dF/dR = n{n + 1) p 2( - n +V / '(- R) n+2 , and, 
using the fact R = —6(H+2H 2 ) at the background level, 
write Eos, pjl - 12"?]) as 



dxi 
dln(a) 

dx 2 
dln(a) 



-1 — 3^2 — X3 + X4 + x\ — x\xz] (39) 
X1X3 



x 2 {xi - 2x 3 + 4) + 



■m 



(40) 
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FIG. 1: (Color online) A phase portrait of the model. The tra- 
jectories in the 3 dimensional space spanned by x\,x^,xz are 
projected to the x\ — X2 and x\ — xz planes respectively. The 
critical points P2,Pz are shown. Only 2 different trajectories 
with different initial conditions are depicted for clearness and 
n = — O.f is adopted. It is seen that x\ oscillates around 0. 
Note that these are just for illustrations and do not represent 
realistic cosmologies. 



dx 3 



eHn(a) 
dln(a) 



2x 3 (2 - x a ) 
x 4 (xi - 2x 3 ). 



XlX 3 



in 



(41) 
(42) 



The dynamics of this system has been thoroughly in- 
vestigated in [3(J, here we just list the main results for 
the model f(R) = R + ^ n+1 '>/(-R) n for the purpose of 
completeness. Firstly, it is easy to obtain 



-R 



n+l 



X 3 + X 2 

x 3 - nx 2 



(43) 



for the use in the numerical calculation. Secondly, the 
critical points and their properties of the system are listed 
in Table [T] where we have defined 



Weft 



ZFH 2 
2H 



= 1 - X\ — X 2 — X 3 — .T 4 , 



1. 



(44) 
(45) 



It is easy to see from Table U that the points Pi, P 2 , P 3 
correspond to a radiation-dominated era, a matter- 
dominated era and a de Sitter era, respectively, while 
P4, P5, Pg cannot lead to any of these three eras and are 
excluded from further consideration. In Table ITTI we sum- 
marise the stability and acceleration properties of the 
critical points other than P4, P5, Pg (limited to the cases 
in which — 1 < n < 0). From this table we see that, for 
general — 1 < n < 0, we can expect Pi,P 2 ,p5 to give 
a viable cosmology provided that the initial conditions 
are chosen appropriately. In particular, around P 2 , the 
evolutionary trajectory of the model is a damped oscil- 
lation, and it is finally attracted to the stable de Sitter 
point P3 . In Figure Q] we have plotted a phase portrait 
to illustrate this evolution: we see that evolutions with 
rather different initial conditions fall towards the vicinity 
of P 2 and finally pass it, being attracted to the de Sitter 
point P3 along nearly the same trajectory (with this tra- 
jectory itself depending on the value of n). Depending on 
the initial conditions, the region of the oscillation could 
be very tiny; but the numerical simulation shows that, 
the oscillation of x\ around is a general result, at least 
for some period of the evolution. 

One can recognise here the connection with our Ein- 
stein frame analysis. Note that ip oc InP, so ip oc F/F 
and an oscillation of x\ around corresponds to an os- 
cillation of F around 0, which in turn represents an os- 
cillation of (p around (this is just the condition for a 
ip oscillation). Furthermore, the final attraction of the 
trajectories towards P3 is also consistent with our con- 
clusion above, that the universe will finally evolve into a 
de Sitter stage. 

In order to go further, and solve the background equa- 
tions numerically, we need the detailed initial conditions. 
One possible way to obtain these is to consider the evo- 
lution in the radiation-dominated era. For example, it 
is possible (in Einstein frame) that the field is frozen 
in this era [4l| while settling in the vicinity of its ef- 
fective potential minimum and then starting to oscillate 
as soon as the fractional matter energy density becomes 
significant. In our calculation, we use a trial and er- 
ror method to find the initial conditions at some speci- 
fied early time which reproduce the observational value 
K/O m o/3i?o — 0-3 an d K/9 ra do/3-ffo — 10~ 4 , where sub- 
script '0' denotes the present-day value (note that we 
are not using f2 m o ~ 0.3, f2 ra do — because of the 

different definitions Eqs. (|36l 22} from conventional cos- 
mology). The results for some specific choices of n are 
given in Figure O Although the difference is not large in 
the plots, we can see that the f(R) dark energy starts to 
be significant earlier for larger |n|'s. This is as expected 
because the ratio of correction term in f(R) to — R is 
given by [/J 2 /{— P)]" +1 which becomes of order 1 earlier 
if n + 1 is closer to 0. On the other hand, evolutions of 
the fractional energy densities are essentially the same as 
that for ACDM at early times when the f(R) corrections 
are negligible. For reference, we also show the case of 
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TABLE I: The critical points and their properties. 



Point 



(a?i, x 2 ,x 3 , Xj) 



Pi 

Pz 

P 3 
Pi 

Ps 
Pe 
Ft 
P« 
P 9 



(0,0,0,1) 

(0,-1,1,0) 

(0,-1,2,0) 
(-1,0,0,0) 
(1,0,0,0) 



4(n + l) 



2(n + l) 
t73 , 



2(n + l) _5n^+8n+2 

77 ' 77^ 



3(77 + 1) 



4n+3 4n+3 

277^ ' 277 



2(77 + 2) 



477 + 5 



,0J 

n(4n + 5) 



2n+l ' (2n+l)(n+l) 5 (2n + l)(n+l) ' ' 

(-4,5,0,0) 




1 

2 



877 2 + 13n+3 






1 






5ti 2 +8ti+2 
TV 



1 
■i 
1 

3 

377 + 4 

3n 

77 + 1 
77 

6ti 2 +7t7-1 
' 3(2ti+1)(t7+1) 
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FIG. 2: (Color online) The fractional energy densities f2 m (green lines), f2 ra d (blue lines) and Qde = 1 — fim — f2rad (red lines) 
as functions of redshift z. We plot the cases for n = —0.01, —0.1, —0.2 and —0.5 respectively. Note that in general Q m Q > 0.3 
because we require Kp m o/3HQ = 0.3 instead and because Fo < 1. 



ACDM (n = 0) in Figure H Secondly, from Eq. (36]) we get 

Next, we outline the way to obtain all other back- x F H 2 

ground quantities from what has already been calculated. — pjj2 a ' 

Firstly, from Eq. (|4"3"1) we now have X40 

so that 

F = (1 + n)— ^ . (46) H 2 = E^a-^H 2 (47) 

X3 — HX2 Fxa 
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TABLE II: The stability and acceleration properties of the critical points (— 1 < n < 0). 



Point 



Eigenvalues 



Stability 



Acceleration ? 



Pi 
Pi 

P 3 
Pe 
Pi 
P« 



1,4,4,-1 



1, 3 + 



3ra 
(1 + m) 2 ' 

3,-4,- 



3 , \/25+if 
2 2 



1^0+ 



1, 



^/81n 2 + 132n + 36 
2n 



-1,3(1 + 1' 



3(n+l)±^/256n 4 +864n 3 + 1025n 2 +498n+81 
4n(n+l) 

4n + 5 2(n + 2) 2(5n 2 +8n + 2) 8n 2 +13n + 3 

n+1 ' 2n + l ' (2n+l)(n+l) ' (2n + l)(n + l) 



Stable if 



Saddle 
Saddle 

Stable 
Saddle 

I < n < ^~ 13 , saddle otherwise 



Stable if n > 13 , saddle otherwise 

— 16 ' 



No 
No 

Yes 
No 
No 

-Kn<-\ h 



"As n — > — 1 this gives a matter era, but one of the eigenvalues 
becomes oo so that the system cannot stay around the point for 
a long time. Also in this case f(R) = R + a/(—R) n oc R, mean- 
ing that the modified f(R) gravity is merely GR with a different 
gravitational constant. 

'Strong phantom era with ui fj < —7.6; unstable. 




Ig(1+z) lg(1+z) 



FIG. 3: (Color online) The same as FIG. 2, but for the case 
of n = (ACDM). 



From Eqs. 



it follows that 



F 
1? 
H 



—xiFH, 
{x 3 - 2)H 2 . 



(48) 
(49) 
(50) 



Finally, Eqs. (|25l [26J) give 

F = (l + 2x 1 +3x 2 +x 3 -x i )FH 2 . (51) 

Consequently, provided Hq is known, we can compute 
any other background quantity. 

The value of F plays an important role in the f(R) 
cosmology and so we discuss it separately here. Dur- 
ing the matter-dominated era, the trajectory oscillates 
around the point P2, at which F = 1 exactly. The rela- 
tion xi + x 3 — holds approximately so that F remains 
close to 1. At the de Sitter attractor P3, we have X2 = — 1 



FIG. 4: (Color online) F as a function of redshift. The 
red, green, blue and magenta curves are the cases for n — 
-0.01, -0.1, -0.2 and -0.5 respectively). 



and X3 — 2, so that F — (2n + 2)/(n + 2) (corresponding 
to the potential minimum at tpo in the Einstein frame). 
Since our universe has not yet reached P3, we expect that 
(2n + 2)/(n + 2) < Fq < 1 at present, which is confirmed 
in Figure [4] where we plot the time evolution of F for 
the same choices of n as above. 



Thus far we have calculated the background quantities 
xi, X2, X3, X4, at some predefined time grid points (scale 
factor values). Then, in the CAMB code, we can obtain 
the background quantities at arbitrary times (scale factor 
values) simply by interpolating between these grid points. 
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FIG. 5: (Color online) The CMB power spectra of the f{R) 
model. The black, red, green, blue and magenta lines repre- 
sent the cases of n — 0,-0.01,-0.1,-0.2 and —0.5 respec- 
tively. For all the curves we have adopted the parameters 
Ho = 70 km/s/Mpc and Kp m0 /^H§ = 0.3. 



C. The CMB and Matter Power Spectra 



FIG. 6: (Color online) The potentials <f>k (see definition in 
text) as functions of the scale factor a for various k val- 
ues. Only the cases of ACDM (n = 0) (black curves), 
n — —0.2 (red curves) and n = —0.5 (green curves) 
are shown for clearness. The scales are chosen to be 
k = 0.0005,0.001,0.005,0.01 Mpc" 1 as indicated beside the 
curves. For all the curves we have adopted the parameters 
Ho = 70 km/s/Mpc and Kp m o = 0.3. 



In the metric f(R) gravity theory, F and R are de- 
termined by the energy-momentum tensor T a b through a 
dynamic equation, as are the perturbations of them, V a F 
or V a i?. To see this, take the covariant spatial derivative 
of the structural equation and we get 

K (V a p-3V a p) = 3(V a F + 9V a F + FV a 6 + V a V 2 F) 



-RV a F - 4r^aF, 

PR 



(52) 



which should be added to our set of perturbation equa- 
tions. Since this is a second-order differential equation, 
we can recast it into two first-order differential equations, 
which means that we have two more quantities to evolve 
in the CAMB code, whose initial conditions also need to 
be specified. 

Making a harmonic expansion of V a F as 



where Q k a = f V a Q k and Q k are the zero-order eigenval- 
ues of the comoving Laplacian a 2 V 2 so that a 2 V 2 Q k = 
k 2 Q k [32j], it is easy to show that 



VaF 

V a F 



(53) 
(54) 



with the aid of 



V a V = (V„V0" + g0V o ^ 



for any scalar ip (the = means that we obtain this relation 
in the frame where the acceleration A a = 0, see |12|). At 
early times, when F = 1, F = 0, Eq. (fS"2")) reduces to 



^ 3 , 



R 
3~ff2 



1 - 



rf r 



U2 1 



a 2 H 2 



e,(55) 



in which X,X P are defined through the harmonic ex- 
pansions V a p = J2k kXQ k /a and V a p = J2k kX p Q k /a, 
and a star denotes the derivative with respect to In a. 
Recall that in the matter-dominated era \R/H 2 \ ~ 1, 
k(X-3X p )/3H 2 = n Pm A m /3H 2 - A m , where A m is the 
fractional perturbation in non-relativistic matter compo- 
nents (the contribution from photons and massless neu- 
trinos to this term is zero because X p v — Xj jV /3) and 

1 1 



RF R 



n(n + l) (4) 



71+1 



tends to —00 (for —1 < n < 0) as /i 2 /(—R) —> 0. So, early 
in the matter-dominated era, the coefficient in front of e 
becomes very large (remember again that R < 0) and e 
settles to ~ A m /|l — 1/RFr\ within a short time [33l |. 
This means that the actual result will be insensitive to 
the choice of initial conditions, a fact we have checked 
using the numerical code. Thus, for our calculation, we 
can simply set einjtiai = einitiai = 0. In this model, the 
perturbation in F is driven essentially to zero (compared 
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with the matter perturbation) and finally grows as dark 
energy starts to drive the expansion of the universe. It 
is also interesting to note that according to the analysis 
above, when n — 0~ we have Fr — _ and so the pertur- 
bation, e, will stay zero all through cosmic history which 
is consistent with the property of an effective cosmologi- 
cal constant. 

In Figure we plot the CMB power spectrum of the 
model for different choices of n. For all of these plots we 
again adopt Ho — 70 km/s/Mpc and rep m o/3i?Q = 0.3. 
Two effects can immediately be seen to occur in the 
spectrum. Firstly, the locations of the acoustic peaks 
move towards the lower Is as \n\ increases. This is be- 
cause, for larger |n|, the dark-energy term starts to be 
important earlier and, with Hq fixed, the universe has a 
smaller age today, which causes the peaks to shift left- 
wards. This effect is not very significant when n is small. 
Secondly, we see a reduction of power at low is, which 
was also been found and discussed in [33j , although there 
the background expansion is fixed to match the ACDM 
cosmology. This reduction in low-/ power is caused by 
a weaker late-time decay of the large-scale potential <f>k 
(which is the coefficient of the harmonic expansion of £ a b 
as £ a b — —J2kk 2 4>kQab/ a2 an d related to the Newto- 
nian potential \& by W = <f> — nlia 2 /2k 2 for any specified 
fc-mode where II is the anisotropic stress [32]) compared 
with that in ACDM, that leads to a suppression on the 
Integrated Sachs- Wolfe (ISW) effect and thus of the C; at 
low Is. In Figure[6]we give some plots for the evolution of 
the potential at different scales (fc's) in this model, which 
clearly show the slower decay at large scales [46]. Note 
that this is remarkably different from the cases of Palatini 
f(R) gravities, in which the potentials decay much more 
rapidly than that in ACDM [lj, , leading to signifi- 
cant amplifications of the low-/ power. Also, we see that 
this model is potentially useful in reducing the difference 
between the predicted and measured CMB powers at low 
I. 

The linear matter power spectra at z = of the present 
models are shown in Figure [7] From this plot we can see 
that the matter power spectrum has similar behavior to 
those in the case of a general f(R) model with ACDM 
background evolution [331 ] - On small scales the spectra 
arising in the f(R) gravity theories we are considering 
have similar shapes to the case of the ACDM power spec- 
trum, and this fact can be understood roughly, as follows. 
Consider for simplicity the growth of the dark-matter 
fractional density perturbation A c , which is defined by 
p c A c = X c (essentially the S c in CMBFAST), in a uni- 
verse filled with cold dark matter and photons. After 
some manipulations of our set of perturbation equations 
it is easy to show that in metric f(R) theories: 



A 



F' 

H ~F 



F 



-nXa 



k 2 e 



F 



He' 



3 H" 



F" 



H' 



2F H 
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FIG. 7: (Color online) The matter power spectra of the f(R) 
model. The black, red, green, blue and magenta lines repre- 
sent the cases of n = 0, —0.01, —0.1, —0.2 and —0.5 respec- 
tively. In all these plots we have adopted Ho = 70 km/s/Mpc 
and Kpmo/SHo = 0.3. 



H 



EL 

2F 



(56) 



where a prime represents the derivative with the confor- 
mal time r (adr — dt) and H = a' /a. At the same time, 
Eq. (|52"1) can now be rewritten as 



k(X - 3X p )a 2 = 3(e" + 2He + kF'Z + k 2 e) 



~6{H' + H 2 ) 1 



RFr) 



e (57) 



where Z is the harmonic expansion coefficient of V a # 

(V a 6< = Efc fc2z Q"/« 2 )- In thc CDM framc ( in which 
v c = 0) we have [321 ] 



A'. 



-kZ. 



(58) 



A similar analysis to that below Eq. ([55)1 then shows that, 
for small scales which are characterised by k 2 jH 2 ^> 1, 
the term 3fc 2 e dominates the right-hand side of Eq. ([57]) 
provided that \Fr\ is not too close to 0, which is the case 
for the later stages in our f(R) cosmology. Thus, we have 
now 



n(X-3X p )a 2 w 3(kF'Z + k 2 e), 



(59) 



which we have also verified numerically. 

Substituting Eqs. {HI [59]) into Eq. ([S6j). we get (ne- 
glecting the contribution from relativistic energy compo- 
nents): 



X + HK-— n Pc A c a 2 



0. 



(60) 



Interestingly, we have obtained a scale-independent be- 
havior for the CDM density perturbation growth on small 
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scales, similar to that of standard GR, where it is given 
by 

A" + HA' C — ^np c A c a 2 » 0. (61) 

This explains why on small scales the shape of the f(R) 
matter power spectrum is like that of ACDM (see [47| 
for a similar example) . It also indicates that the growths 
of small-scale perturbations are governed by an effective 
gravitational constant given by K e ff = 4re/3.F > k. Notice 
that one should not simply try to put F = 1 into Eq. (|6"0|) 
to get the ACDM limit, since under this condition the 
last term in Eq. (|57|) always dominates over 3k 2 e and can 
never be neglected. The correct ACDM limit is recovered 
by setting F' = 0,F = 1,6(H' + H 2 )(l - F/RF R )e = 
-k(X - 3X p )a 2 , (according to Eq. and all other 

terms relating e and e' to in Eq. (fB"6")) . which just leads 
to Eq. (|6"Tj) , and shows that the early- matter-era growth 
of perturbations is well described by GR. 

We can now make a comparison between the matter 
power spectra in metric and Palatini f(R) cosmologies. 
In the later, it has been shown that any small deviation 
from GR will cause the small-scale spectra to differ from 
the GR one in obscrvationally unacceptable ways, either 
blowin g up or oscillating and being prevented from grow- 
ing (Ta,ll4l[l5|. This is because in the Palatini f(R) grav- 
ity F and T satisfy an algebraic equation F = F(T) such 
that those terms involving V 2 F can always be rewritten 
as cx V 2 T ~ V 2 p and consequently a term cx k 2 A c ap- 
pears in the growth equation which makes the growths on 
small scales strongly scale-dependent. In the metric f(R) 
gravity here, in contrast, e is determined by the perturba- 
tions in T through a differential equation Eq. ([57)1 . As k 
increases, the value of the perturbation in F(e) decreases, 
so that k 2 e does not exceed 0(kXci 2 ) and becomes k in- 
dependent, so its effect cannot be as exotic as that arising 
in Palatini f(R) theories. 

Before leaving this analysis, we want to briefly dis- 
cuss the cosmological viability of the present model. A 
more rigorous analysis involves carefully searching the 
multidimensional parameter space (as in [3) and is be- 
yond the scope of the present work. Firstly, we have 
seen that the background evolutions allowed by the model 
are rather close to the ACDM paradigm with particular 
choices of n, and so could be consistent with the SNIa ob- 
servations. The confrontation between predictions and 
the data on the linear power spectrum is a little more 
complicated. For the CMB spectrum, Figure O indicates 
that it is similar to the ACDM prediction at higher I's 
and can also reduce the quadrupole power and bring it 
closer to the measured values [33[ (We have also checked 
that this low-? power reduction is a general feature of the 
model and is insensitive to the values of Hq and p m o). 
However, because of the limitation from cosmic variance, 
the influences on the likelihood analysis are expected to 
be small. For the matter spectrum, at small scales the 
power is significantly higher than for the ACDM case, 
yet the fact that it has similar shape to the latter indi- 



FIG. 8: (Color online) The dependence of the matter power 
spectra on Ho (upper panel) and Kp m o/3Ho (lower panel). 
Upper panel: the black, red, green solid curves represent Ho = 
65,68,70 km/s/Mpc respectively (Kp m o/3i?o = 0.3). Lower 
panel: the black, red, green curves represent fcp m o/3iifo = 
0.25,0.27,0.30 respectively (H = 70 km/s/Mpc). Dashed 
curves are the spectra for the ACDM model, n = —0.5 is 
adopted for illustration. We see that a smaller Kp m o/3-ffo can 
bring the shape of the spectrum closer to the ACDM one: in 
this case data on large scale structure cannot place stringent 
constraints on the model. 



cates the possibility that a different choice of galaxy bias 
might be able to make the model's predictions consistent 
with, for example, the Sloan Digital Sky Survey (SDSS) 
observations. However, we notice that at wavenumbcr 
k ~ 0.05 /iMpc^ 1 deviations from the ACDM shape start 
to be significant, which might potentially be in conflict 
with the data; this could still be alleviated by reducing 
the quantity Kp m o/3ffg, as can be seen in Figure [H 

We conclude that constraints on this model from large- 
scale structure could be much weaker than those on the 
Palatini f(R) models Il5| . This said, it is still in- 
teresting to obtain quantitative joint constraints from 
background cosmology, CMB, and matter power spec- 
tra on the model, and examine whether it can be made 
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consistent with other cosmological observations (such as 
considered in (33| ) in the future. 

IV. DISCUSSION AND CONCLUSIONS 

To summarise: in this article we have considered the 
cosmology of a subclass of metric f(R) gravity theories, 
that are characterised by f(R) = R + fj. 2( . n+1 ^ / (- R) n , 
and have discussed both the flat background Friedmann 
universe and its inhomogeneous perturbations. For the 
background evolution, we address the problem in both 
the Jordan frame and the Einstein frame and find the cor- 
respondences between them. Generally, the evolution is 
attracted towards a saddle point in the phase space which 
has the characteristics of a matter-dominated era. If the 
initial conditions are chosen appropriately, the universe 
can stay in the vicinity of that point for a sufficiently 
long period of CDM-dominated evolution to occur. Fi- 
nally, it always passes this saddle point and evolves to 
a stable de Sitter point. The cosmic expansion begins 
to accelerate during this transition period and a cosmic 
history that is consistent with observations is possible in 
this model, as expected from (30| . In the Einstein frame, 
the model reduces to one with a scalar field conformally 
coupled to (non-relativistic) matter, and it shares some 
of the properties of a chameleon field cosmologies. That 
cosmological model results in a matter-dominated era fol- 
lowed by an accelerating era can also be seen from the 
analysis in this frame. The ACDM limit of the model is 
also discussed. 

We derive the covariant and gauge invariant perturba- 
tion equations for f(R) gravity theories in general, which 
we believe to be rather convenient for numerical calcula- 
tions, and applied them to the present model through a 
modified CAMB code. The linear spectra are obtained 
and discussed. The CMB power spectrum displays reduc- 
tions in the low-Z power, which arises from a weakening 
of the ISW effect because of insufficient late-time decay 
of the large-scale potentials, as shown in the plots of the 
evolution. For the matter power spectrum, we find that 
on small scales the matter density growth is nearly scale- 
independent, making the shape of the spectrum at large 
fc's similar to that of the ACDM spectrum. Compar- 
isons with the CMB and matter power spectra in Palatini 
f{R) gravity theories are then made, which account for 
the dramatic differences between these two approaches 
to modify GR, despite the fact that they could have the 
same gravitational action. 

Possible comparisons with different observational data 



sets are also briefly discussed. It is found that neither the 
data on background evolution, CMB spectrum, nor those 
on the matter power spectrum can be used to exclude the 
model. Their joint constraints on the model, however, 
could be complicated and involve making a numerical 
search of the parameter space, which is beyond the scope 
of the present work. However, we can see that constraints 
on this metric f(R) model could be much weaker than 
those on the Palatini f(R) theories, because for the later 

(1) the CMB power is largely amplified at low Vs and 

(2) the matter power spectrum at small scales is strongly 
scale dependent, both conflicting with observations on 
CMB anisotropies and large scale structure. 

The form of f(R) adopted here is one of the few which 
could produce a viable model for the entire cosmic his- 
tory. In this model the modifications to GR take effect 
only at late times. Therefore, it is interesting to look for 
other models in which the modifications are significant at 
different times, along the lines of refs. [l5|, [3J]. In [HI, 
a model with f(R) = R + Xi exp(i?/A2) is considered in 
the Palatini approach and its effects on the linear power 
spectra are also discussed. In (34|, this same model is 
investigated within the metric approach; however, it ac- 
tually gives a <f>- matter-dominated era (0MDE), [24L|48|. 
rather than the standard matter-dominated era, and is 
not cosmologically viable. Thus, it may still be meaning- 
ful to look for viable early-time f(R) cosmologies in the 
metric formalism. Another interesting issue to explore is 
whether the f(R) gravity models might be more adapted 
to hot dark matter. These topics will be investigated 
elsewhere. 
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